Variational Mean Field approach to the Double Exchange Model 
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It has been recently shown that the double exchange Hamiltonian, with weak antiferromagnetic 
interactions, has a richer variety of first and second order transitions than previously anticipated, 
and that such transitions are consistent with the magnetic properties of manganites. Here we 
present a thorough discussion of the variational Mean Field approach that leads to the these results. 
We also show that the effect of the Berry phase turns out to be crucial to produce first order 
Paramagnetic-Ferromagnetic transitions near half filling with transition temperatures compatible 
with the experimental situation. The computation relies on two crucial facts: the use of a Mean 
Field ansatz that retains the complexity of a system of electrons with off-diagonal disorder, not fully 
taken into account by the Mean Field techniques, and the small but significant antiferromagnetic 
superexchange interaction between the localized spins. 
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I. INTRODUCTION 

Doped manganites show many unusual features, the 
most striking being the colossal magnetoresistance 
(CMR) in the ferromagnetic (FM) phase In addi- 

tion, the manganites have a rich phase diagram as func- 
tion of band filling, temperature and chemical composi- 
tion. The broad features of these phase diagrams can 
be understood in terms of the double exchange model 
(DEM) ||,|, although Jahn-Teller deformations D and 
orbital degeneracy may also play a role 0|. A remark- 
able property of these compounds is the existence of in- 
homogeneities in the spin and charge distributions in a 
large range of dopings, compositions and temperatures 
[p[-pi)[. In fact, for materials displaying the largest CMR 
effects, the size of the phase-separated domains is so large 
(~ 0.5/im Q), that the electrostatic stability of the ma- 
terial should be addressed by theorists. At band fillings 
where CMR effects are present, x ~ 0.2 — 0.5, these com- 
pounds can be broadly classified into those with a high 
Curie temperature and a metallic paramagnetic (PM) 
phase, and those with lower Curie temperatures and an 
insulating magnetic phase pd]-p^. 

The double exchange mechanism was introduced by 
Zenner Q through the following Kondo lattice type 
model 



(1) 



where t and Jh are, respectively, Cg electron's hopping 
and Hund's coupling between the Cg and the localized 
t2g electrons responsible for the core spin S. 

When Jh is larger than the width of the conduc- 
tion band, the model can be reduced to the double ex- 
change model with weak inter-atomic antiferromagnetic 
(AFM) interactions. Early investigations showed a 
rich phase diagram, with AFM, canted and FM phases, 



depending on doping and the strength of the AFM cou- 
plings. More recent studies have shown that the compe- 
tition between the double exchange and the AFM cou- 
plings leads to phase separation into AFM and FM re- 
gions, suppressing the existence of canted phases [|5|-^. 
In addition, the double exchange mechanism alone in- 
duces a change in the order of the FM transition, which 
becomes of first order, and leads to phase separation, at 
low doping Note, however, that a systematic study 
of the nature of the transition at finite temperature was 
not addressed until recently despite its obvious rel- 
evance to the experiments. In fact, in Ref. it was 
shown that a small AFM uniform superexchange inter- 
action between the localized t2g spins is crucial to un- 
derstand some of the more relevant features of the phase 
diagram of the manganites. In particular a first-order 
phase transition is found between the PM and FM phases 
in the range x ^ 0.2 — 0.5 is found. This transition does 
not involve a significant change in electronic density, so 
that domain formation is not suppressed by electrostatic 
effects. Therefore, we find a phase-separation of a rather 
different type of the previously discussed, not driven by a 
charge instability, but by a magnetic instability. In addi- 
tion to this phase transition, we recover those previously 
discussed. 

In this work we give a detailed exposition of the new 
mean-field technique |20| and emphasis is made on the 
importance of the Berry phase for the existence of first or- 
der phase transitions near half filling. We have been able 
of achieving a more complete description of the phase di- 
agram than in previous work, because we have taken full 
profit of a very particular feature of the DEM, namely, 
fermions are bilinearly coupled to classical degrees of 
freedom (the Mn spins). This allows to trace-out the 
fermions, thus obtaining a non local effective Hamilto- 
nian for the spins, that can be explicitly written down in 
terms of the density of states of the fermionic hopping 
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matrix. What we propose to do is to calculate exactly 
the effective spin-Hamiltonian, for a given (disordered in 
general) configuration of the spins, using the so called 
moments-method |2l|| complemented with an standard 
truncation procedure . This technique can be directly 
used for other models, like for instance models of classi- 
cal spins and lattice vibrations coupled to fermions with- 
out direct interactions or also in contexts different 
from the manganites physics like the pyrochlores or dou- 
ble perovskytes. Once electrons are traced-out, we study 
the spin thermodynamics using the variational version of 
the Weiss Mean-Field method ^ . 

The main difference of our approach with previous 
work is in that we use the exact spin-Hamiltonian, while 
an approximated form was used up to now. For instance, 
the effective crystal approximation, which amounts to 
consider that electrons move on a perfect crystal, with a 
magnet ically reduced hopping, was employed in p4| (see 
section VIII for a comparison with our method) . A more 
accurate estimate of the density of states, well-known 
from the physics of disordered-systems and which be- 
comes exact on the infinite dimensions limit p9| ], is the 
CPA. Notice that for non self-interacting electrons the 
Dynamical Mean-Field Approximation(DMFA) is 
also equivalent to a CPA calculation on a given mean- 
field, which is determined via self-consistency equations 
(both the mean-field and the CPA density of states are 
obtained self-consistently) . Although those approaches 
are reasonable, and provide useful information, they are 
approximated in two different ways, which is undesider- 
able because two different effects are entangled. First, 
even for only-spins models of magnetism (like the Ising 
Model) where the exact spin Hamiltonian can be evalu- 
ated easily, the Mean-Field approximation neglects the 
spatial correlations of the statistical fluctuations of the 
order parameter (see e.g. |2^). But, in addition, for an 
electron-mediated magnetic interaction, the evaluation 
of the effective spin Hamiltonian is only accurate in the 
limit of infinite dimensions. With our approach, the cor- 
relations on the magnetic fiuctuations are to be blamed 
for all the differences between our results and the real be- 



havior of the model. On the other hand, in section VIII 



we show how the failure of the effective-crystal approxi- 
mation on finding the first-order phase transition at half- 
filling is due to the inaccuracy on the calculation of the 
density of states. Moreover, we are able to study directly 
in three dimensions some rather subtle details, like the 
non-negligible effects of keeping the Berry phase on the 
DEM Hamiltonian. 

The structure of the paper is as follows. In Section II 
we introduce the DEM and our notations. In Section HI 
we present our Mean-Field approximation. The very non- 
trivial part of the work, the computation of the Density 
of States, is explained in Section IV. It requires numer- 
ical simulations that can be performed on large lattices 
with a high accuracy. They are described in Section V. 
The effects of the Berry phase are analyzed in Section VI. 
Section VII is devoted to the study of the influence of the 



Berry phase in the phase diagram of the DEM. The com- 
parison of the Mean-Field approach studied in this work 
with the de Gennes's jl^ and with the DMFA is carried 
out in Section VIII. The conclusions are summarized in 
Section IX. 



II. MODEL 



We study a cubic lattice with one orbital per site. At 
each site there is also a classical spin. The coupling be- 
tween the conduction electron and this spin is assumed 
to be inflnite, so that the electronic state with spin an- 
tiparallel to the core spin can be neglected. Finally, we 
include an AFM coupling between nearest neighbor core 
spins. We neglect the degeneracy of the conduction band. 
Thus, we cannot analyze effects related to orbital order- 
ing, which can be important in the highly doped regime, 
X > 0.5 |]7|,p6[ (sec however Ref. [^). We also neglect 
the coupling to the lattice. We focus on the role of the 
magnetic interactions only. As mentioned below, mag- 
netic couplings suffice to describe a number of discon- 
tinuous transitions in the regime where CMR effects are 
observed. These transitions modify substantially the cou- 
pling between the conduction electrons and the magnetic 
excitations. Thus, they offer a simple explanation for 
the anomalous transport properties of these compounds. 
Couplings to additional modes, like optical or acoustical 
phonons, will enhance further the tendency towards first 
order phase transitions. We consider that a detailed un- 
derstanding of the role of the magnetic interactions is re- 
quired before adding more complexity to the model. Note 
that there is indeed evidence that, in some compounds, 
the coupling to acoustical phonons |^ or to Jahn- Teller 
distortions is large. 

The Hamiltonian of the DEM is 



JafS Si ■ Sj 



(2) 



where S" = 3/2 is the value of the spin of the core, Mn'^+, 
and S stands for a unit vector oriented parallel to the 
core spin, which we assume to be classical. In the fol- 
lowing, we will use Jaf = JafS'^. Calculations show 
that the quantum nature of the core spins does not in- 
duce significant effects In one of the earliest studies 
of this model , the superexchange coupling was cho- 
sen FM between spins lying on the same z = constant 
plane, and AFM between spins located on neighboring 
planes. This is a reasonable starting point for the study 
of Lai_a;Ca2;Mn03 if a; < 0.16, where A-type antiferro- 
magnetism is found. For larger doping, 0.16 < x < 0.5, 
which is our main focus, the magnetism is uniform and 
there is no a priori reason for favoring a particular direc- 
tion. 

The function 
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cos — COS 
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sm - 



(3) 



stands for the overlap of two spin 1/2 spinors oriented 
along the directions defined by Si and Sj, whose polar 
and azimuthal angles are denoted by and ip, respec- 
tively. It defines a hopping matrix T, whose matrix ele- 
ments are Tij — T{Si, Sj). The hopping function can be 
written as 



'^{Si, Sj 



icos-^exp(i(?!)y), 



(4) 



where dij is the relative angle between Si and Sj, and 
(pij is the so called Berry phase. It is sometimes assumed 
that the Berry phase can be set to zero without essential 
loss. It is therefore interesting to study the model that 
ignores the Berry phase, the hopping matrix being 



t:"-'^ ^17^,1= cos % 



Si • Si 



(5) 



In the following sections we will analyze both models, 
with Berry phase (hopping matrix T) and without Berry 
phase (hopping matrix 7"™°'^). 



III. MEAN-FIELD APPROXIMATION 

Our approach to the problem follows the variational 
formulation of the Mean-Field approximation, described 
for instance in [Q. We start by writing the Grand 
Canonical partition function for the DEM: 



--GC 



(6) 



where ^ is the electronic chemical potential, M — c^Ci 
is the electron number operator, T is the temperature 
and we use units in which the Boltzmann constant is one. 
The trace, taken over the electron Fock space, defines an 
effective Hamiltonian for the spins, 

exp [-'W^{S)/T] = Tr(P°='^) exp [-(W - ^iN)lT] , (7) 

that can be computed in terms of the eigenvalues. En, of 
the hopping matrix, T: 

H^^^iS) = Jaf ^ Si ■ Sj - 

T^log{l + cxp[-(i?„(5)-M)/T]} . (8) 

n 

Introducing the Density of States (DOS) of T: 

g{E;S)^-J2S[E-En{S)], (9) 



where V is the volume of the lattice, the effective Hamil- 
tonian can be written as 

H^^{S) — ^ Jaf Si ■ Sj ~ 



TV / dEgiE;S) log 



(10) 



The Grand Canonical partition function becomes an in- 
tegral in spin configuration space: 



2gc - J [dS] cxp[-Heff(5)/r] 



(11) 



Thermodynamics follows from Eq. ( [1 iD as usual. The 
free energy, and the electron density, x, are given by: 



T 
V 



log Zgc 



'^^JdEimS))- 



. exp [{E - m)/T] 



(12) 



(13) 



where (• ■ ■) stands for spectation value over equilibrium 
spin configurations. 

The variational Mean-Field approach consist on com- 
paring the actual system with a set of simpler reference 
models, whose Hamiltonians, Tih, depend on external pa- 
rameters, hi. For simplicity, we choose the model: 



J2hi-Si 



The variational method is based on the inequality 



(14) 



(15) 



where J-h is the free energy of the system with Hamil- 
tonian (|l^), and the expectation values (• • ■)h are calcu- 
lated with the Hamiltonian Tih. The inequality ( p^ fol- 
lows easily from the concavity of the exponential function 
p3[ . The best approximation to the actual free energy 
with the ansatz of Eq. (|l4[) is 



.F = min{j-^ + (W^ff)^-(7iO^} • 



(16) 



Since, for technical reasons that will become clear in 
the following, it is not possible to work with one field hi 
per site, we must select some subsets that contain only a 
few independent parameters (see Section VII) . The choice 
is of paramount importance since it is an ansatz that will 
artificially restrict the behavior of the system. We have 
chosen the following four families |^ of fields, depending 
on a parameter, h: 



hi^h, 

K = {-ly-h, 
h, = i-ir+y^h, 

hi = (-l)^'+2''+^-/i, 



(17) 
(18) 
(19) 
(20) 
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which correspond, respectively, to FM, A-AFM, C-AFM, 
and G-AFM orderings. There is an order parameter 
(magnetization) associated to each of these orderings. 
We win denote them by Afp, Ma, Mq, and Mq, respec- 
tively. As a shorthand, they will be denoted generically 
by M ■ The order parameter is related to the correspond- 
ing h by 



M = 



1 



tanh h 



1 

h 



(21) 



where h = \h\/T. Thus, the free energy can be written 
in terms of A4 instead of h and Eq. ( [iq ) implies that it 
must be minimized with respect to A4. The free energy 
has three contributions: the fermion free energy (FFE), 
the superexchange energy and the entropy of the spins: 



TiM) - TfUM) + NJafM^ - TSh{M) 



(22) 



where N is , respectively, 3,-3, 1 and —1 for FM, G- 
AFM, A-AFM and C-AFM orderings. 

The entropy of the spins can be easily computed in 
terms of the Mean-Field: 



Sh{M) = {Th-{nh)H)lT 

= log [sinhh{M)/h{M)] - h{M)M , 



(23) 



but the FFE, 



^Fcr(A^) =-T j dE {g{E;S)) ^\og 



1 + c-(E-,.)/T 



(24) 



must be estimated numerically. 

The non trivial part of the computation is the aver- 
age of the DOS {g{E; S))h- The key is that it can be 
computed by numerical simulations on large lattices with 
high accuracy, ought to two basic facts: 1) the Mean 
Field Hamiltonian (^4|) describes uncorrelated spins and 
therefore equilibrium spin configurations can be easily 
generated on very large lattices, and 2) the DOS is a self- 
averaging quantity. This last point means that the mean 
value of the DOS can be obtained on a large lattice by av- 
eraging it over a small set of equilibrium configurations. 
Once the DOS is computed, the integral of Eq. ( p^ ) can 
be performed numerically to get the FFE as a function 
of A^. 



Now, it is easy to show that the eigenvalues of the hop- 
ping matrix verify —6t < En < 6t. A probability distri- 
bution of compact support can be reconstructed from its 
moments using the techniques of Stieljes |3l|| . In prac- 
tice, we only know the first p moments, but the method 
of Stieljes allows us to find a good approximation to the 
distribution if p is large enough. 

To compute the averaged DOS we follow four steps: 

i) Generate spin configurations, {5'}, according to 
the Mean-Field Boltzmann weight, exp(— Ti^/T")- 
This can be achieved very efficiently with a heat 
bath algorithm, since all the spins are decoupled in 
the Mean-Field Hamiltonian. In this way, one ob- 
tains spin configurations in perfect thermal equi- 
librium with the Boltzmann-weight given by the 
Mean-Field Hamiltonian. 

ii) For each {S}, one would calculate the moments of 
the DOS, using Eq. (|25|), and then apply the tech- 
niques of reconstruction of Stieljes pl[ . However, 
this is impractical since, although the matrix T is 
sparse, the trace in Eq. ( |25| ) would require to repeat 
the process V times, and we would end-up with an 
algorithm of order . We use instead an stochastic 
estimator. First, we extract a normalized random 
vector \v), with components 



(26) 



where the ai are random numbers extracted with 
uniform probability between — 1 and 1 . Let us now 
call \n) to the eigenvector of eigenvalue i?„ of the 
matrix T. It is easy to check that (the overline 
stands for the average on the random numbers ai) 



ViV 



{n\v){v\ra 



' - V ' 



V 



(27) 
(28) 



Then we introduce a w-dependent density of states: 

V 

g{E;v;S) = Y,\{M^)\'SiE-En). (29) 



IV. COMPUTATION OF THE AVERAGED DOS 

The DOS can be accurately computed for any given 
spin configuration with the technique that we describe 
in the following ||2^. From its definition, Eq. (^, the 
DOS is a probability distribution in the variable E, whose 
moments are 



MS) = / dEgiE;S)E'' 



1 
V 



-TrT'^ 



(25) 



From Eq. (|28|), it follows immediately that 



giE;v;S)^g{E;S). 



(30) 



From {v\v) = 1, and from Eq. (p9|), we see that 
g{E; v; S) is a perfectly reasonable distribution 
function, whose moments are 



J dE E''g{E-v-S) = {v\T>^\ 



(31) 
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Numerically, the algorithm is of order k xV, since, 
as mentioned, T is sparse and only 0{V) opera- 
tions are required to multiply v by T. This method 
allows to compute a large number of moments on 
large lattices. However, notice that the actual cal- 
culation is not performed this way (round-off errors 
would grow enormously with the power of T), but 
as explained in the appendix. 

iii) Reconstruct g{E;v;S) from the moments, by the 
method of Stieljes. The DOS is obtained in a dis- 
crete but very large number of energies, E. The 
cost of refining this set of energies is negligible. 
Hence, the integral over E that gives the FFE, 
Eq. (24), can be approximated numerically with 



s(i') ; 

0.25 — 



high accuracy. 

iv) Average g{E; v; S) over the spin configurations S 
and over the random vectors \v). In practice we 
only use a random vector per spin configuration: 
it is useless to obtain an enormous accuracy on the 
density of states for a particular spin configuration, 
that should be spin-averaged, anyway. Since the 
errors due to the fluctuations of the spins and the 
fluctuations of the \v) are statistically uncorrelated, 
both of them average out simultaneously |Q . It is 
crucial that the DOS is a self-averaging quantity, 
what means that its fluctuations are suppressed as 
1/VV. Hence, its average over a few equilibrium 
configurations is enough to estimate it with high 
accuracy. 

This program can be carried out successfully on lat- 
tices as large as 64'^ and even 96"^, computing a large 
enough number of moments, 50 or 75. As we shall see, 
this suffices to achieve an excellent accuracy in the aver- 
aged DOS and in the FFE. 

The whole process is repeated for several values of the 
Mean Fields h within each family. The FFE computed in 
a discrete set of magnetizations is extended to a contin- 
uous function of A4 in the interval [0, 1] through a poly- 
nomial fit, as we will show in next section. The magnetic 
phase diagram of the model will come out easily then. 



EXTRACTING THE FERMIONIC FREE 
ENERGY 



Let us discuss in this section the numerical results for 
the averaged DOS and for the FFE. We carried out the 
program designed in the previous section for 20 values of 
h within each family of Mean Fields, chosen in such a way 
that the corresponding magnetizations, M. , are uniformly 
distributed in [0, 1]. The expectation value {g{E\ S) )h is 
estimated by averaging over 50 equilibrium (with respect 
to the Mean Field distribution) spin configurations. This 
is enough to have the statistical errors under control on 
a 64^ lattice. 
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FIG. 1. Averaged DOS, reconstructed with 50 moments, 
versus E, for four values of the mean-field corresponding to 
My = 1 (solid line), Mp = 0.5, (dashed line), paramagnetic 
(dotted line), and Mq — 0.5, (dot-dashed line). 

Figure |^ displays the averaged DOS, computed on a 
64'^ lattice for four values of h, corresponding to FM, 
(Mp — 0.5 and Alp = 1), paramagnetic, and G-AFM 
{Mq = 0.5) phases. They were reconstructed with its 
50 first moments Note that the DOS is even in 

E, as required by the the particle-hole symmetry, and 
that it becomes narrower in going from the FM to the 
G-AFM phase, as expected. The width of the density 
of states in the PM phase is 8t, two thirds of the width 
corresponding to the perfect ferromagnetic, in agreement 
with the results of full diagonalization Q| . 

From the averaged DOS it is straightforward to 
compute the FFE by performing the integral entering 
Eq. (^) numerically. In this way, we computed J^fci- 
for the chosen values of A4. Given that the free en- 
ergy can be shifted by a term independent of A4, we 
use J^Fcr(A^) - ^For(O) instead of J^fcAM) ||. 

To have an analytic expression for J-por, we fit the data 
with a polynomial of order sixth in A4 , with coefficients 
that depend on T and /i: 



42(M/) = A'^^'M'f + A'i'Mf + A\'' M'j. (32) 

The index / denotes the type of ordering: / ~ F, A, C, 
or G. 

Figure ^ shows the FFE at T = and /i = 0, which 
correspond to half filling, x — 1/2, as a function of Mf 
or Mq. The points are the result of the numerical com- 
putation and the lines are the best fits of the form (|3^). 
The high quality of the fits is remarkable. Note that the 
fermions favor FM order. The coefficients A'"p{T,ii) of 
Eq. (^), which will play a major role in the exploration 
of the phase diagram, are displayed in Figure |^ in the 
cases / = F, G, for T = 0, as a function of ^. We always 
found Af\T,^Ji) < and A^^\t,^jl) > 0, in agreement 
with the FM nature of the spin interaction induced by 
the double exchange mechanism. 
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FIG. 2. FFE at T = and ^ = (x = 1/2) versus Mp 
(squares, solid line) or Mg (circles, dotted line). The points 
are the results of the simulation and the lines are the three 
parameter fit. 
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FIG. 3. Coefficients A^P , Af\ Af\ Af \ Af\ and Af' 
of the fits (^^ versus ^, At T — 0. Each figure contains four 
lines (that sometimes cannot be distinguished) corresponding 
to the four different computations mentioned in the text, three 
with 50 and one with 75 moments. 

Let us estimate the errors of our numerical approach. 
We have three sources of errors: 

a) Finite size of the lattice. 

b) Statistics, arising from the numerical simulation. 

c) Truncation of the infinite sequence of moments. 

Finite size errors have been estimated comparing the 
results on a 64^ and a 96'^ lattice. They turn out to be 
negligible, as expected given the sizes of the lattices. To 
estimate the statistical errors we performed three differ- 
ent simulations using the 50 lowest order moments. One 
more simulation, this time with the 75 first moments, 
was done in order to study the systematic error associ- 
ated to the truncation of the sequence of moments. As 



an example. Figure |4| displays the averaged DOS in the 
FM phase (Mp = 0.5) extracted from two different sim- 
ulations, one with 50 and the other with 75 moments. 
We see only tiny differences, which can hardly be ap- 
preciated on the scale of the figure. This small error 
propagates to the FFE, which is shown in Figure ^ for 
T = and ji — 0. There are four sets of points plot- 
ted, corresponding to the four mentioned simulations. 
Again, the differences cannot be appreciated. The er- 
rors in the computed FFE give rise to uncertainties in 
the coefficients A's of Eq. (^2[), which can be appreci- 
ated in Fig. ||. The largest errors appear at /z = (half- 
filling). We have also checked that fits to a polynomial of 

eight-order do not change the values of A^p significantly. 




FIG. 4. Averaged DOS vs. energy for Mp = 0.5 on a 
64^ lattice from a simulation using 50 moments (solid line) 
and another one with 75 moments (dashed). The curves can 
hardly be distinguished on this scale. 
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FIG. 5. Fermion free energy vs. Mf from four simulations 
on a 64^ lattice. Three of them, carried out to estimate the 
statistical errors, used 50 moments to reconstruct the DOS. 
The other one, aimed to estimate the systematic errors due to 
the truncation of the sequence of moments, took into account 
75 moments. The errors turn out to be so small that cannot 
be appreciated on the scale of the figure. 



To summarize, we have checked that the numerical un- 
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certainties inherent to our numerical approach are well 
under control and therefore all the conclusions are robust 
in this sense. The discrepancies between our analysis and 
the true behavior of the system, if they are important, 
must be attributed only to the Mean-Field ansatz. 



VI. EFFECT OF THE BERRY PHASE 

Sometimes it is stated that one can ignore the Berry 
phase in the DEM without loosing anything important. 
To investigate this point, let us repeat our analysis of the 
DEM setting the Berry phase to zero. All we have to do 
is to compute the DOS of the hopping matrix 7^™°*^ of 
Eq. (|^). The rest is identical to what we have discussed 
in the previous sections. 



g 
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FIG. 6. Averaged DOS in the paramagnetic phase 
(Mp = 0, upper panel) and in the ferromagnetic phase 
(Mf = 0.5, lower panel) with (solid) and without (dashed) 
Berry phase. 

Figure | displays the averaged DOS at Mp = (PM 
phase) and M-p = 0.5 for hopping with and without Berry 
phase. At first sight, the differences, although noticeable 
(they are much bigger than the errors, cf. Fig. I), do 
not seem very important. However, it happens that the 
results are very sensitive to small modifications of the 
DOS. We shall see indeed that the presence or absence of 
the Berry phase is crucial for some features of the phase 
diagram. Thus, the analysis of errors of the previous 
section turns out to be extremely important to give a 



meaning to our results. Notice that the effect of the Berry 
phase is stronger in the disordered phase, as expected. 
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FIG. 7. Fermion free energy at T = and fj. — as a func- 
tion of Mf with (solid) and without (dashed) Berry phase. 
Notice that, in both cases, the free energy has been taken 
vanishing at the origin by convention. 

Figure shows the effects of the Berry phase on the 
FFE at T = and ^ = 0. These effects modify the coef- 
ficients A^i^^ entering J-per, which can be seen in Fig. 
These coefficients are very sensitive to the modifications 
of the DOS induced by the Berry phase. Of especial 
relevance is A'"^\ which, as we shall see in the next sec- 
tion, governs the possibility of having first order PM-FM 
transitions. In particular, in the vicinity of /i = this 
coefficient is negative. Notice however that without the 
Berry phase A^^"* is negative around ^ — in a, smaller 
region than with Berry phase, and it is closer to zero. 
This fact induces important differences in the nature of 
the phase transitions of the model, as we shall see in the 
next section. 



VII. PHASE DIAGRAM 

The equilibrium states are determined by the absolute 
minima of the free energy, Eq. (|2^), respect to the or- 
der parameters, A4. The minima determine the phases 
and the phase boundaries. Given that we know !F as 
a function of Ai, the problem of determining the equi- 
librium states is reduced to numerical minimization of a 
function of a single variable. This is indeed the way we 
proceed. It is however illuminating to get some insight 
by a semi-analytic treatment of the problem. As we have 
seen, to a very good approximation, the FFE is a poly- 
nomial of sixth degree in Ai. The entropy ( p3| ) can also 
be expanded in powers of A4 around = 0: 

+ + + ...). (33) 
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Hence, we find the Landau expansion of the free energy 
in powers of the order parameter: 



T{M) = C2M^ + CiM^ + ceM^ + 



(34) 



iii) Degeneracy (:^(0) = T{Mo)): 

C2 + CiM'o + ceMl 



(41) 
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FIG. 8. Coefficients of the fit ( p2| ) of the fermion free energy 
at T = as a function of fj. with (solid) and without (dashed) 
Berry phase. 



The coefficients of the expansion are 

C2 = ^T + NJaf + 4^^(r,M), 
ca^^T + A(^)(T,,), 



(35) 
(36) 
(37) 



where N was defined right after Eq. (|2S 

The free energy has the symmetry Ai — > . At high 
T, the entropic term dominates and the minimum of J-' is 
at = 0. As the temperature decreases, the internal en- 
ergy becomes more important and the absolute minimum 
of J- can be located at 0. The phase transition will 

be continuous when the absolute minimum at the origin 
changes to a maximum, i.e.: 



C2 







+ A^Jaf + A^2\Tc,^i) = (38) 



At a first order transition three minima, A4 = and 
A4 — ±A4q 0, are degenerate. Three conditions must 
hold: 



i) Minimum at = 0: 

C2 > 0. 

ii) Minimum at A^q- 

2c2 + AC4MI + GceMt = . 



(39) 



(40) 



T/t 
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FIG. 9. Phase diagram of the DEM in the plane {x,T/t), 
for several values of JAp/t. The left (right) part corresponds 
to the model with (without) Berry phase. Solid (dashed) lines 
represent first (second) order transitions and the zones with 
stripes are phase separation regions. The onset for first order 
PM-FM transition is at Jaf ~ 0.06 in the model with Berry 
phase, while such transitions do not appear if the Berry phase 
is neglected. 



The solution of these three equations is: 



Mt 



Cq — '-.4 
C2 > 



= — 2C2/C4 
C4/(4C2) 



C4 < 0; cg > 



(42) 
(43) 
(44) 



Eq. (|4^) gives the spontaneous magnetization; Eq. (|42 
determines the critical temperature of the first order 
transition as a function of and Jaf; Eq. (|^ sets nec- 
essary conditions (real Alo) for a first order transition to 
happen. We see that to have a first order PM-FM tran- 
sition we must have A^P (Tc, /i) < 0. Fig. || shows that in 
particular this is possible around half filling. 

The boundaries between first and second order lines are 
tricritical points. They are determined by the conditions: 



C2 = 
C4 = 



P — 



3 

Jt — 



2N 



N 



20 



T, = -=-Ai'\T„fi) 



(45) 
(46) 
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With these ingredients, we are able to discuss the 
phase diagram of the DEM. It has been shown in |2^] 
that the phase diagram of double exchange systems is 
richer than previously anticipated and differs substan- 
tially from that of more conventional itinerant ferro- 
magnets. Moreover, it is consistent with the magnetic 
properties of manganites p|- [T3| , ^6| -^ (see section IX). 
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00 
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C-AFM - 













FM 







0,0 0,1 0,2 0,3 0,4 0,5 



FIG. 10. Phase diagram of the DEM without Berry phase 
in the plane {x, Jaf) at T = 0. 

We shall not repeat here the analysis of the phase dia- 
gram of the DEM carried out in |^ . Let us concentrate 
on the effects of the Berry phase. Figure || displays the 
phase diagram in the plane (x, T/t), for several values of 
Jaf- The left part corresponds to the model that includes 
the Berry phase and on the right the Berry phase is ne- 
glected. For Jaf < 0.06 both phase diagrams are very 
similar. The transition temperature is slightly higher 
(13%) at half filling in the model that ignores the Berry 
phase. At low filling, both phase diagrams are almost 
identical. However, for Jaf > 0.06 important differences 
arise. A first order PM-FM transition develops around 
half filling if the Berry phase is taken into account. The 
onset for such behavior is Jaf = 0.06. In contrast, the 
PM-FM transition around half filling remains continuous 
for any value of Jaf if the Berry phase is neglected. The 
explanation is simple: the coefficient A\ , although neg- 
ative, has a too small absolute value to drive a first order 
transition. This negative coefficient would become ef- 
fective if the critical temperature were much lower, what 
can be achieved by increasing the AF superexchange cou- 
pling. But in this case the competition between FM and 
AF is so strong that the transition at half filling takes 
place between PM and A-AFM phases, and it is second 
order. First order PM-FM transitions only appear if the 
Berry phase is properly taken into account. 

At T = 0, the phase diagrams are similar in both cases, 
and we only display that of the model without Berry 
phase, in Fig. |l^. The discussion of applies to this 
case without any modification. 



Let us end this section with the analysis of the phase 
transitions at finite applied magnetic field, B. The first 
order PM-FM transition around half filling survives un- 
der an applied magnetic field. In this case, the order 
parameter, M-p, is non-zero in both phases, but suffers 
a jump on a line in the plane (B,T). The line ends at 
a critical point, {B*,T*), which has a certain magneti- 
zation Mp. The critical field can be measured and is of 
interest p9|p0[ |. Let us compute it. The free energy in 
the presence of a magnetic field, B, is: 

T{Mf) ^ C2M| + C4M| + ceM§ - BMp . (47) 

The magnetic field shifts the three degenerate minima 
of the zero-field PM-FM first order transition and lifts 
the degeneracy. By tuning (increasing) the temperature 
it is possible to get two degenerate minima again, and 
a first order transition takes place. In this way, we get 
a transition line in the [B, T) plane. Increasing B, the 
two degenerate minima become closer. At the critical 
field, B* , both minima coalesce at some point Mp, and 
the transition disappears. When this happens, the three 
first derivatives of J- respect M-p vanish, and the fourth 
is positive. These conditions read: 

T'{M^) = 0^ B* = 2c2Mf + 4c4M^3 ^ QceMf^ (48) 

T"iM^) = ^ 2c2 + 12c4M^2 ^ 30c6M;4 = (49) 

T"'iM^) = ^ 24c4M; + 120c6M;3 ^ Q (50) 

J^(™)(A/;) > ^ 24c4 + SGOcqM;'^ > 0. (51) 

These equations determine B* , Tc and Mp as a function 
of /i (or x) and Jaf- The critical temperature Tc varies 
very little from its value at i? = 0. Fig. |ri| displays B*, 
in units of 10"'* t, versus x, for Jaf = 0.08 In physical 
units, using t « 0.166eV, B* varies from 0.6 Tesla at 
X = 0.33 to 2.2 Tesla at a; = 0.5. Recent measurements 
in Lao. 6 Y0.07 Cao.33 MnOs gave a critical field of 1.5 Tesla 




0,25 0,30 0,35 0,40 0,45 0,50 

X 

FIG. 11. Critical magnetic field (in units of 10"* t) versus 
X at Jaf = 0.08 1. 
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VIII. COMPARISON WITH OTHER 
CALCULATIONS 

A. Rigid band mean field approximation. 

The main conclusion of the variational Mean-Field 
technique applied to the DEM is the prediction of a first 
order PM-FM transition at half-filling and its vicinity 
for JAp/t S [0.06, 0.1]. This is in sharp contrast with the 
widely used Mean-Field approach devised by de Gennes 
in 1960 |l4j, which predicts a second order PM-FM at 
half- filling for any value of Jaf- Let us see briefly what 
are the differences between these two approaches that 
yield different qualitative behavior. 

The difficulty in the Mean-Field approach to the DEM 
lies in calculating the contribution of the fermions to the 
Mean-Field free energy. In the variational method dis- 
cussed here, we compute it exactly through a numerical 
simulation. As we have already mentioned, the only ap- 
proximation is the Mean-Field ansatz for the Boltzmann 
weights of the spin configurations. On the other hand, 
de Gennes suggested that the fermion free energy might 
be well approximated by the free energy of an assem- 
bly of fermions propagating on a crystal with an homo- 
geneous hopping parameter given by the average of the 
spin-dependent hopping parameter over the Mean-Field 
spin configurations. The de Gennes' method neglects 
the influence of the Berry phase. In this approach, the 
electronic DOS depends on the spin configuration only 
through the hopping parameter (in this case, without 
Berry phase): g{E; S) = g{E;T'^°'^{Si ■ Sj)). In mathe- 
matical terms, de Gennes's approximation is carried out 
through the following substitution: 

{g{E; |r(S. • S,)\))^ ^ giE; {\T{S, ■ S,)|)^) 

= go{E;To) (52) 

where go{E; Tq) is the DOS of free fermions with hopping 

%ih) = {\T{S,-S,)\)h 



1 



6 



18 



^i/2H/^)S(2;-l)(2Z + 3)' 



J 



l+l/2\ 



-ih) 



(53) 



and Ji,{z) is the Bessel function. 

Since the hopping is homogeneous, the fermion free 
energy is known analytically. At T = and half-filling 
{fi — 0) it is: 



dEg{E-%)E 



-ToWo 



Wo 



dEgoiE;l)E 



(54) 
(55) 



All the dependence in the magnetization is contained in 
To{h). The expansion in powers of the magnetization fol- 
lows straightforwardly from Eqs. (B3) and (21). It yields: 



175 



875 



(56) 



The coefficient of Mp in JFpor is positive. Hence, the PM- 
FM phase transition at half-filling can only be continu- 
ous. We have also checked that this remains true when 
we keep the contribution of all powers of Mp to J-pcr- 

The fermions in de Gennes's approach propagate only 
on perfect crystals. In the truly variational Mean-Field 
presented in this work, the fermions propagate on the 
disordered spin background generated by the Mean-Field 
h. This appears as an important ingredient that leads 
the predictions closer to the phenomenology, as we have 
shown in Ref. fS. 



B. Dynamical Mean Field Approximation. 

This method allows for an improvement on the treat- 
ment of the electronic contribution to the self energy. In 
the PM phase, the density of states is proportional to 
that in the fully ferromagnetic phase, like in de Gennes 
treatment. The only difference is that the constant of 
proportionality is l/\/2 and not 2/3. Below Tc, the den- 
sity of states is calculated self consistently, through a self 
energy which can be written as: 



and: 



^{E;S^,)^{\T{S,-S,)\'g{E;S,)}s^ 



E-^{E- S, 



(57) 



(58) 



Finally, the average (■ • ■)« is carried out defining a prob- 
ability distribution, V{Sj), which depends self consis- 
tently on the free energy associated with a site with 
magnetization Sj immersed in the lattice described by 

Our approach is similar to the dynamical mean field 
approximation, but differs from it in two aspects: 

i) The electronic density of states is calculated in a 
cubic lattice, instead of using the semielliptical DOS valid 
in the Bethe lattice with infinite coordination. 

ii) We use a variational ansatz for V{Sj), instead of 
determining it fully self consistently. 

Point i) allows us to consider effects of the lattice geom- 
etry, and the influence of the Berry's phase, as discussed 
above. At zero temperature, where both approaches be- 
come exact for their respective lattices, we find phases 
which can only be defined in a 3D cubic lattice. 

If the transition is continuous, the distribution V{Sj) 
can be expanded on the deviation from the isotropic one, 
'P{Sj) = const. , in the PM phase. The ansatz that 
we use has the correct behavior sufficiently close to Tc, 
so that both approaches will predict the same value of 
Tc, for a given lattice. One must be more careful in 
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the study of discontinuous transitions. Our ansatz in- 
troduces an approximation in the ordered phase (which 
disappears at T — 0). However, near a first order tran- 
sition we do not expect divergent critical fluctuations, 
so that our approach should give qualitative, and prob- 
ably semiquantitative correct results as compared to the 
DMFA, in lattices where the latter is exact. 



C. Hierarchy of approximations. 



We are tempted to design a hierarchy of approxima- 
tions, ordered according to the cocfhcient ^4 of the Alp 
term in the Landau expansion of the fermion free energy, 
as follows: 

(F) 

1) de Gennes's approximation: ^4 > and the PM- 
FM transition is second order. 

2) Exact variational computation without Berry 



phase: ^4^^ < but \A''^'\ too small to produce 
first order PM-FM transitions, see Eq. (|36|). 

3) Exact variational computation with Berry phase: 



(F), 



A\ ' < and IA4 | large enough to produce first 
order PM-FM transitions. 



IX. CONCLUSIONS 

We have presented a detailed analysis of the variational 
Mean Field technique. This method can be useful in any 
situation where non self-interacting fermions are coupled 
to classical continuous degrees of freedom. Within this 
method, the fermionic contribution to the free energy is 
calculated exactly, and, later on, the variational Mean 
Field method is applied to the classical degrees of free- 
dom. As an example, we have chosen the Double Ex- 
change Model, both with and without Berry phase. The 
phase diagram has been obtained in both situations. 

We have shown that the Berry phase is crucial in order 
to get first order PM-FM phase transitions around half 
filling. Such transitions are second order if the topolog- 
ical effects associated to the Berry phase are neglected. 
Thus, the dimensionality of the lattice plays a very im- 
portant role in the structure of the phase diagram. 

Some earlier Mean-Field computations approx- 
imate the fermion free energy by that of an assembly of 
fermions propagating on a perfect crystal with an ho- 
mogeneous hopping parameter averaged over the spin 
configurations. They yield second order FM-PM tran- 
sitions in the vicinity of half-filling. The propagation 
of the fermions in the disordered spin background gen- 
erated by the Mean-Field is another crucial ingredient 
to get discontinuous PM-FM transitions at half-filling. 
More modern approaches, as the DMFA ||2^ cannot deal 
with three dimensional effects as the Berry phase either. 



As shown in |Q , the variational Mean-Field described 
in the present work leads to results that are consistent 
with the phenomenology of the magnetic properties of 
the manganites La.i-x (Sr, Ca)^ MnOs, in the range 0.3 < 
X < 0.5, in particular with the fact that for materials with 
a high transition temperature, the PM-FM transition is 
continuous while for those with low Tc is not. Moreover, 
the order of magnitude of our estimate of the critical 
field for which histeretic effects disappear agrees with the 
experimental findings in Lao.eo Yq.q? Cao.33 MnOs 
Also the phase diagram obtained by substitution of a 
trivalent rare earth for another one with smaller ionic ra- 
dius (i.e. compositional changes that do not modify the 
doping level) is in remarkable agreement with our results. 

Of course, the DEM itself can also be highly improved. 
For instance, one should include the orbital degeneracy, 
which is known to play an important role, and other el- 
ements like phonons and Jahn- Teller distortions. The 
variational Mean-Field approach can be applied with the 
same techniques presented in this work whenever the 
bosonic fields that interact with the electrons can be 
treated as classical. 
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APPENDIX A: THE METHOD OF MOMENTS 

In this appendix we include, for completeness, some 
details on the method of moments For a com- 

plete mathematical background we refer to Ref. [^ . 
The method of moments allows to obtain some statis- 
tical properties of large matrices (as the density of states 
or the dynamical structure factors, in general quantities 
depending on two-legged Green functions) , without actu- 
ally diagonalizing the matrices. Regarding the density of 
states, once one recognizes that it is a probability func- 
tion whose moments can be obtained by iteratively mul- 
tiplying by T the initial random vector \v), it is clear 
that the classical Stieljes techniques |31| can be used. 
Ought to the fact that the matrix T is sparse, and using 
the random vector trick, Eq. (^T]), the moments can be 
calculated with order V operations. Since the spectrum 
of the matrix T lies between —6t and 6t for any spin 
configuration, the Stieljes method is guaranteed to con- 
verge. The procedure is as follows: one first introduce 
the resolvent 



R(z) 



-6t 



E' 



(Al) 
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that has a cut along the spectrum of T, with discontinu- 
ity 

2TTg{E) = Im Urn [R{E - ie) - R{E + ie)] . (A2) 

The resolvent can be obtained from the orthogonal 
polynomials of the g{E), with the monic normalization: 



starting with 



P„(S) = E^ + Cn-iE" 



dE g{E) Pr,{E) P„,{E) . 



(A3) 



(A4) 



with Pq = 1, P_i = 0, and n, m = 0, 1, . . . The polyno- 
mials verify the following recursion relation: 

Pn+l{E) = {E-an)Pn{E) - b,, Pn-l{E) , (A5) 

with the coefficients a„ and &„ given by: 



dEg{E)EPl{E) 



dEg{E)P^{E) 



bn = 



r_l^ dEg{E)Pl{E) 
dEg{E)P?^_,{E) 



(A6) 
(A7) 



The coefficient bo is arbitrary and is conventionally set- 
tled to one. 

The resolvent has a representation in terms of a con- 
tinued fraction as follows: 



Riz) 



1 



z — ao — 



bi 



(A8) 



If one truncates the continued fraction, the resolvent 
would be approximated by a ra tion ale function, which 
does not have a cut and use Eq. ( |A^ ) is impossible. For- 
tunately, when, as in this case, the density of states does 
not have gaps, the coefficients q„ and 6„ tends fastly to 
their asymptotic values a and b ||3l| . Thus, one can end 
the continued fraction [ p2[ with a truncation factor T{z), 
that verifies 



T{z) 



z ^ a — T{z) 



(A9) 



Since the previous equation is quadratic in T{z), we find 
that T[z) has a branch-cut between a — 2\/b and a + 2^/b, 
which are the limits of the spectrum. 

One should not use the moments of the g{E) calculated 
with Eq. ( |3l| ) to obtain the orthogonal polynomials (and 
hence the {a„,6„}), since this is an extremely unstable 
numerical procedure. It is better to use the recurrence 
relation: 

Pn+i{T)\v) = (T - a„)F„(T)|z;) - 6„P„_i (r)|i;) , 

(AlO) 



P-,{T)\v) = 



From this, one immediately gets 

_ {Pr,{T)v\TPr,{T)v) 



bn = 



{P„{T)v\P„iT)v) ' 

{Pr,iT)v\Pr,iT)v) 
{Pn-l{T)v\Pn-l{T)v) 



(All) 

(A12) 
(A13) 



In this way, one generates the A^*'' orthogonal polyno- 
mial of the matrix (times v) recursively, at the price of 
N multiplications per T. The cost of this procedure is 
always of order V operations. For each random-vector, 
one first extracts the density of states through Eq. ( [A2| ) , 
which is subsequently averaged over the different \v) and 
spin realizations. 

Let us finally point out that the above recursion rela- 
tion is virtually identical to the Lanzcos method (the only 
difference lies on the normalizations). It should therefore 
not be pursued for a large number of orthogonal poly- 
nomials, without reorthogonalization. For the relative 
modest number of coefficients calculated in this work, 
this has not been needed. 
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